---
title: "Public Perceptions of ‘Woke’ Corporate Political
Advocacy"
author: "Wayde Z.C. Marsh & Jordan Carr Peterson"
date: "Last Edited: 02/07/26"
output:
  html_document:
    toc: true
    code_folding: hide
  latex_engine: xelatex
  number_sections: yes
mainfont: Helvetica
editor_options: 
  chunk_output_type: inline
---

# Setting up the Data {.tabset}

## Loading packages and data

```{r Packages, message = FALSE, warning = FALSE, results ='hide'}

##### Packages #####

library(descr)
library(dplyr)
library(tidyverse)
library(ggplot2)
library(ggpubr)
library(ggeffects)
library(data.table)
library(haven)
library(pwrss)
library(lmtest)
library(qualtRics)
library(modelsummary)
library(nnet)
library(sandwich)
options("modelsummary_format_numeric_latex" = "mathmode")

# set working directory

# setwd(dirname(rstudioapi::getSourceEditorContext()$path)) 
# knitr::opts_knit$set(root.dir = '~/')

```

## Clean Study 1 Data

```{r Clean Study 1 Data}

dat1 <- read_survey('CPT2_Study1_raw.csv')

# convert dataframe to base r dataframe object
dat1 <- as.data.frame(dat1)

# prune respondents who failed attention checks
dat1 <- subset(dat1, attn1 == 'Yellow' & attn2 == 'Strongly agree')

dat1 <- dat1[, c(1:2, 21:102, 109:113)] # remove unnecessary variables

dat1 <- dat1 %>%
  mutate(chevron_exp = recode(chevron_exp, 'chevron_exp1|chevron_timer' = 'Control',
                           'chevron_exp2|chevron_timer' = 'JudicialExpertise',
                           'chevron_exp3|chevron_timer' = 'Regulation'),
   hatespeech_exp = recode(hatespeech_exp, 'hatespeech_exp1|hatespeech_timer' = 'Control',
                           'hatespeech_exp1|hatespeech_timer' = 'AntiImmigrant',
                           'hatespeech_exp1|hatespeech_timer' = 'AntiLGBT'),
         cpt_exp = recode(cpt_exp, 'cpt_exp1|cpt_timer' = 'Control',
                           'cpt_exp2|cpt_timer' = 'Fiscal',
                           'cpt_exp3|cpt_timer' = 'Woke',
                           'cpt_exp4|cpt_timer' = 'WokeFiscal'),
   elect_concerns_exp = recode(elect_concerns_exp, 'elect_concerns_exp1' = 'Control',
                           'elect_concerns_exp2' = 'Treatment1',
                           'elect_concerns_exp3' = 'Treatment2',
                    'elect_concerns_exp1|elect_concerns_exp2|elect_concerns_exp3' = 'NA'),
   sexassault_exp = recode(sexassault_exp, 
                           'sexual_assault_exp1|sexual_timer' = 'WomanControl',
                           'sexual_assault_exp2|sexual_timer' = 'ManControl',
                           'sexual_assault_exp3|sexual_timer' = 'WomanTreatment',
                           'sexual_assault_exp4|sexual_timer' = 'ManTreatment'))

dat1$elect_concerns_exp <- ifelse(dat1$elect_concerns_exp == 'NA', NA, dat1$elect_concerns_exp)

dat1 <- dat1 %>%
  mutate_at(c('cpt_dv1','cpt_dv2'),
            funs(dplyr::recode(., 'Strongly disagree' = 0, 
                           'Disagree' = 0.25,
                           'Neither agree nor disagree' = 0.5,
                           'Agree' = 0.75,
                           'Strongly agree' = 1, .default = NaN)))

dat1 <- dat1 %>%
  mutate(cpt_dv4 = recode(cpt_dv4, 
                          'Punish corporations that engage in this behavior' = 0,
                          'Do nothing regarding corporate behavior' = 0.5,
                          'Reward corporations that engage in this behavior' = 1))
dat1$cpt_dv3 <- dat1$cpt_dv3 / 100

dat1$pid <- trimws(dat1$pid)
dat1$partyid <- 
  ifelse(dat1$pid == 'Democrat' & dat1$dem1 == 'Strong Democrat', 0,
  ifelse(dat1$pid == 'Democrat' & dat1$dem1 == 'Not a very strong Democrat', 1,
  ifelse((dat1$pid == 'Independent' & dat1$ind1 == 'Democratic party') |
         (dat1$pid == 'Something else' & dat1$ind1 == 'Democratic party'), 2,
  ifelse((dat1$pid == 'Independent' & dat1$ind1 == 'Neither party') |
         (dat1$pid == 'Something else' & dat1$ind1 == 'Neither party'), 3,
  ifelse((dat1$pid == 'Independent' & dat1$ind1 == 'Republican party') |
         (dat1$pid == 'Something else' & dat1$ind1 == 'Republican party'), 4,
  ifelse(dat1$pid == 'Republican' & dat1$rep1 == 'Not a very strong Republican', 5,
  ifelse(dat1$pid == 'Republican' & dat1$rep1 == 'Strong Republican', 6, NA)))))))
dat1$pid2 <- ifelse(dat1$partyid < 3, 'Democrat',
             ifelse(dat1$partyid > 3, 'Republican', NA))

# write.csv(dat1, '~/Dropbox/Projects/CPT/Data/CPT2_Study1.csv')
# write.csv(dat1, '/Users/wmarsh5/Dropbox/Wayde/Research/CES 2024/CES2024/Study1.csv')

```

## Clean Study II Data

```{r Clean Study II Data}

dat2 <- read_sav('~/Dropbox/Wayde/Research/CES 2024/UTK/CCES24_UTK_OUTPUT.sav')

# convert dataframe to base r dataframe object
dat2 <- as.data.frame(dat2)

dat2 <- dat2[, c(3, 7, 8, 10, 11, 13, 28, 196, 207, 258, 364, 373:383)]

dat2$age <- 2024 - dat2$birthyr

dat2$sex <- ifelse(dat2$gender4 == 2, 1,
            ifelse(dat2$gender4 == 1, 0, NA))
dat2$religiosity_pew <- ifelse(dat2$pew_churatd == 6, 0,
                        ifelse(dat2$pew_churatd == 5, 1,
                        ifelse(dat2$pew_churatd == 4, 2,
                        ifelse(dat2$pew_churatd == 3, 3,
                        ifelse(dat2$pew_churatd == 2, 4,
                        ifelse(dat2$pew_churatd == 1, 5, NA))))))
dat2$edu <- dat2$educ - 1
dat2$Frame <- ifelse(dat2$UTK311 == 1, 'Control',
                ifelse(dat2$UTK311 == 2, 'Fiscal',
                ifelse(dat2$UTK311 == 3, 'Woke', 'WokeFiscal')))
dat2$cpt_dv4 <- ifelse(dat2$UTK315 == 1, 0,
                ifelse(dat2$UTK315 == 3, 0.5,
                ifelse(dat2$UTK315 == 2, 1, NA)))

dat2$cpt_dv5 <- ifelse(dat2$UTK316_1 == 1, 'Increase taxes',
                ifelse(dat2$UTK316_2 == 1, 'Offer tax breaks/incentives',
                ifelse(dat2$UTK316_3 == 1, 'Publicly applaud the corporations',
                ifelse(dat2$UTK316_4 == 1, 'Publicly criticize the corporations',
                ifelse(dat2$UTK316_6 == 1, 'State legislatures should do nothing', 'Something else')))))

dat2$cpt_dv5_text <- ifelse(dat2$cpt_dv5 == 'Something else', dat2$UTK316_t,
                            dat2$cpt_dv5)

dat2$pid2 <- ifelse(dat2$pid7 < 4, 'Democrat',
             ifelse(dat2$pid7 > 4 & dat2$pid7 < 8, 'Republican', NA))

dat2 <- dat2[, c(1, 31, 7, 8, 23:26, 5, 6, 27, 12:14, 28, 29, 30)]

names(dat2) <- c('teamweight', 'pid2', 'pid7', 'ideo5', 'age', 'sex', 'religiosity_pew', 
                 'edu', 'race', 'hispanic', 'Frame', 
                 'cpt_dv1', 'cpt_dv2', 'cpt_dv3', 'cpt_dv4', 'cpt_dv5', 
                 'cpt_dv5_text')

dat2$cpt_dv1 <- (dat2$cpt_dv1 - 1) / 4
dat2$cpt_dv2 <- (dat2$cpt_dv2 - 1) / 4
dat2$cpt_dv3 <- dat2$cpt_dv3 / 100

# write.csv(dat2, '~/Dropbox/Projects/CPT/Data/CPT2_Study2.csv')
# write.csv(dat2, '/Users/wmarsh5/Dropbox/Wayde/Research/CES 2024/CES2024/Study2.csv')

```

# Load Data

```{r Load Data}

dat1 <- read.csv('CPT2_Study1.csv')
dat2 <- read.csv('CPT2_Study2.csv')

dat1$dv1 <- scale(dat1$cpt_dv1)
dat1$dv2 <- scale(dat1$cpt_dv2)
dat1$dv3 <- scale(dat1$cpt_dv3)
dat1$dv4 <- scale(dat1$cpt_dv4)

dat1 <- dat1 %>%
  mutate(sex = recode(gender, 'Female' = 1,
                           'Male' = 0,
                           'Non-binary / third gender' = 2),
         religiosity_pew = recode(religiosity_pew,
                            'Never' = 0,
                            'Seldom' = 1,
                            'A few times a year' = 2,
                            'Once or twice a month' = 3,
                            'Once a week' = 4,
                            'More than once a week' = 5,
                            "Don't know/prefer not to say" = NaN),
         edu = recode(educ, 
                      'Did not graduate from high school' = 0,
                      'High school graduate (including GED)' = 1,
                      'Some college, but no degree' = 2, 
                      "Associate's Degree (2-year college)" = 3,
                      "Bachelor's Degree" = 4,
                      "Postgraduate degree" = 5))
dat1$sex <- scale(dat1$sex)
dat1$age1 <- scale(dat1$age)
dat1$edu <- scale(dat1$edu)
dat1$religosity_pew <- scale(dat1$religiosity_pew)
dat1$Frame <- dat1$cpt_exp

dat2$sex <- scale(dat2$sex)
dat2$age1 <- scale(dat2$age)
dat2$edu <- scale(dat2$edu)
dat2$ideo <- scale(dat2$ideo5)
dat2$religosity_pew <- scale(dat2$religiosity_pew)
dat2$dv1 <- scale(dat2$cpt_dv1)
dat2$dv2 <- scale(dat2$cpt_dv2)
dat2$dv3 <- scale(dat2$cpt_dv3)
dat2$dv4 <- scale(dat2$cpt_dv4)

```

# Analysis

```{r Analysis}

# interactive model

m1 <- lm(dv1 ~ Frame*pid2, data = dat1)

m2 <- lm(dv2 ~ Frame*pid2, data = dat1)

m3 <- lm(dv3 ~ Frame*pid2, data = dat1)

m4 <- lm(dv4 ~ Frame*pid2, data = dat1)

summary(m1); summary(m2); summary(m3); summary(m4)

m1_2 <- lm(dv1 ~ Frame*pid2, weights = teamweight, data = dat2)

m2_2 <- lm(dv2 ~ Frame*pid2, weights = teamweight, data = dat2)

m3_2 <- lm(dv3 ~ Frame*pid2, weights = teamweight, data = dat2)

m4_2 <- lm(dv4 ~ Frame*pid2, weights = teamweight, data = dat2)

summary(m1_2); summary(m2_2); summary(m3_2); summary(m4_2)

# subgroup analysis

m1_d <- lm(dv1 ~ Frame, data = subset(dat1, pid2 == 'Democrat'))
m1_d2 <- lm(dv1 ~ Frame, weights = teamweight, data = subset(dat2, pid2 == 'Democrat'))

m1_r <- lm(dv1 ~ Frame, data = subset(dat1, pid2 == 'Republican'))
m1_r2 <- lm(dv1 ~ Frame, weights = teamweight, data = subset(dat2, pid2 == 'Republican'))


m2_d <- lm(dv2 ~ Frame, data = subset(dat1, pid2 == 'Democrat'))
m2_d2 <- lm(dv2 ~ Frame, weights = teamweight, data = subset(dat2, pid2 == 'Democrat'))

m2_r <- lm(dv2 ~ Frame, data = subset(dat1, pid2 == 'Republican'))
m2_r2 <- lm(dv2 ~ Frame, weights = teamweight, data = subset(dat2, pid2 == 'Republican'))

m3_d <- lm(dv3 ~ Frame, data = subset(dat1, pid2 == 'Democrat'))
m3_d2 <- lm(dv3 ~ Frame, weights = teamweight, data = subset(dat2, pid2 == 'Democrat'))

m3_r <- lm(dv3 ~ Frame, data = subset(dat1, pid2 == 'Republican'))
m3_r2 <- lm(dv3 ~ Frame, weights = teamweight, data = subset(dat2, pid2 == 'Republican'))

m4_d <- lm(dv4 ~ Frame, data = subset(dat1, pid2 == 'Democrat'))
m4_d2 <- lm(dv4 ~ Frame, weights = teamweight, data = subset(dat2, pid2 == 'Democrat'))

m4_r <- lm(dv4 ~ Frame, data = subset(dat1, pid2 == 'Republican'))
m4_r2 <- lm(dv4 ~ Frame, weights = teamweight, data = subset(dat2, pid2 == 'Republican'))


summary(m1_d); summary(m1_r)
summary(m2_d); summary(m2_r)
summary(m3_d); summary(m3_r)
summary(m4_d); summary(m4_r)


# With Age and Edu controls

# interactive model

m1_c <- lm(dv1 ~ Frame*pid2 + age1 + edu, data = dat1)

m2_c <- lm(dv2 ~ Frame*pid2 + age1 + edu, data = dat1)

m3_c <- lm(dv3 ~ Frame*pid2 + age1 + edu, data = dat1)

m4_c <- lm(dv4 ~ Frame*pid2 + age1 + edu, data = dat1)

summary(m1); summary(m2); summary(m3); summary(m4)

m1_2_c <- lm(dv1 ~ Frame*pid2 + age1 + edu + ideo, weights = teamweight, data = dat2)

m2_2_c <- lm(dv2 ~ Frame*pid2 + age1 + edu + ideo, weights = teamweight, data = dat2)

m3_2_c <- lm(dv3 ~ Frame*pid2 + age1 + edu + ideo, weights = teamweight, data = dat2)

m4_2_c <- lm(dv4 ~ Frame*pid2 + age1 + edu + ideo, weights = teamweight, data = dat2)

summary(m1_2); summary(m2_2); summary(m3_2); summary(m4_2)

# subgroup analysis

m1_d_c <- lm(dv1 ~ Frame + age1 + edu, data = subset(dat1, pid2 == 'Democrat'))
m1_d2_c <- lm(dv1 ~ Frame + age1 + edu + ideo, weights = teamweight, data = subset(dat2, pid2 == 'Democrat'))

m1_r_c <- lm(dv1 ~ Frame + age1 + edu, data = subset(dat1, pid2 == 'Republican'))
m1_r2_c <- lm(dv1 ~ Frame + age1 + edu + ideo, weights = teamweight, data = subset(dat2, pid2 == 'Republican'))


m2_d_c <- lm(dv2 ~ Frame + age1 + edu, data = subset(dat1, pid2 == 'Democrat'))
m2_d2_c <- lm(dv2 ~ Frame + age1 + edu + ideo, weights = teamweight, data = subset(dat2, pid2 == 'Democrat'))

m2_r_c <- lm(dv2 ~ Frame + age1 + edu, data = subset(dat1, pid2 == 'Republican'))
m2_r2_c <- lm(dv2 ~ Frame + age1 + edu + ideo, weights = teamweight, data = subset(dat2, pid2 == 'Republican'))

m3_d_c <- lm(dv3 ~ Frame + age1 + edu, data = subset(dat1, pid2 == 'Democrat'))
m3_d2_c <- lm(dv3 ~ Frame + age1 + edu + ideo, weights = teamweight, data = subset(dat2, pid2 == 'Democrat'))

m3_r_c <- lm(dv3 ~ Frame + age1 + edu, data = subset(dat1, pid2 == 'Republican'))
m3_r2_c <- lm(dv3 ~ Frame + age1 + edu + ideo, weights = teamweight, data = subset(dat2, pid2 == 'Republican'))

m4_d_c <- lm(dv4 ~ Frame + age1 + edu, data = subset(dat1, pid2 == 'Democrat'))
m4_d2_c <- lm(dv4 ~ Frame + age1 + edu + ideo, weights = teamweight, data = subset(dat2, pid2 == 'Democrat'))

m4_r_c <- lm(dv4 ~ Frame + age1 + edu, data = subset(dat1, pid2 == 'Republican'))
m4_r2_c <- lm(dv4 ~ Frame + age1 + edu + ideo, weights = teamweight, data = subset(dat2, pid2 == 'Republican'))

```

## Tables

```{r Tables}

# Interaction Models

modelsummary(models = list('CPT--Climate'  = m1,
                           'CFS' = m2,
                           'CPT' = m3,
                           'Punish--Reward'  = m4),
             output = 'latex',
             stars = c('*' = .05),
             gof_omit = 'AIC|BIC|F|RMSE|Log.Lik.',
             escape = FALSE)

modelsummary(models = list('CPT--Climate'  = m1_2,
                           'CFS' = m2_2,
                           'CPT' = m3_2,
                           'Punish--Reward'  = m4_2),
             output = 'latex',
             stars = c('*' = .05),
             gof_omit = 'AIC|BIC|F|RMSE|Log.Lik.',
             escape = FALSE)

modelsummary(models = list('CPT--Climate'  = m1_c,
                           'CFS' = m2_c,
                           'CPT' = m3_c,
                           'Punish--Reward'  = m4_c),
             output = 'latex',
             stars = c('*' = .05),
             gof_omit = 'AIC|BIC|F|RMSE|Log.Lik.',
             escape = FALSE)

modelsummary(models = list('CPT--Climate'  = m1_2_c,
                           'CFS' = m2_2_c,
                           'CPT' = m3_2_c,
                           'Punish--Reward'  = m4_2_c),
             output = 'latex',
             stars = c('*' = .05),
             gof_omit = 'AIC|BIC|F|RMSE|Log.Lik.',
             escape = FALSE)

# Party Subgroup Models

modelsummary(models = list('CPT--Climate'  = m1_d,
                           'CFS' = m2_d,
                           'CPT' = m3_d,
                           'Punish--Reward'  = m4_d,
                           'CPT--Climate'  = m1_r,
                           'CFS' = m2_r,
                           'CPT' = m3_r,
                           'Punish--Reward'  = m4_r),
             output = 'latex',
             stars = c('*' = .05),
             gof_omit = 'AIC|BIC|F|RMSE|Log.Lik.',
             escape = FALSE)

modelsummary(models = list('CPT--Climate'  = m1_d2,
                           'CFS' = m2_d2,
                           'CPT' = m3_d2,
                           'Punish--Reward'  = m4_d2,
                           'CPT--Climate'  = m1_r2,
                           'CFS' = m2_r2,
                           'CPT' = m3_r2,
                           'Punish--Reward'  = m4_r2),
             output = 'latex',
             stars = c('*' = .05),
             gof_omit = 'AIC|BIC|F|RMSE|Log.Lik.',
             escape = FALSE)

modelsummary(models = list('CPT--Climate'  = m1_d_c,
                           'CFS' = m2_d_c,
                           'CPT' = m3_d_c,
                           'Punish--Reward'  = m4_d_c,
                           'CPT--Climate'  = m1_r_c,
                           'CFS' = m2_r_c,
                           'CPT' = m3_r_c,
                           'Punish--Reward'  = m4_r_c),
             output = 'latex',
             stars = c('*' = .05),
             gof_omit = 'AIC|BIC|F|RMSE|Log.Lik.',
             escape = FALSE)

modelsummary(models = list('CPT--Climate'  = m1_d2_c,
                           'CFS' = m2_d2_c,
                           'CPT' = m3_d2_c,
                           'Punish--Reward'  = m4_d2_c,
                           'CPT--Climate'  = m1_r2_c,
                           'CFS' = m2_r2_c,
                           'CPT' = m3_r2_c,
                           'Punish--Reward'  = m4_r2_c),
             output = 'latex',
             stars = c('*' = .05),
             gof_omit = 'AIC|BIC|F|RMSE|Log.Lik.',
             escape = FALSE)

```

## Figures

```{r Figures}

# Study I

coef1 <- c(m1_d$coefficients[2], m2_d$coefficients[2], # dems coefs fiscal
           m3_d$coefficients[2], m4_d$coefficients[2],
           m1_r$coefficients[2], m2_r$coefficients[2], # reps coefs fiscal
           m3_r$coefficients[2], m4_r$coefficients[2],
           m1_d$coefficients[3], m2_d$coefficients[3], # dems coefs woke
           m3_d$coefficients[3], m4_d$coefficients[3],
           m1_r$coefficients[3], m2_r$coefficients[3], # reps coefs woke
           m3_r$coefficients[3], m4_r$coefficients[3],
           m1_d$coefficients[4], m2_d$coefficients[4], # dems coefs wokefiscal
           m3_d$coefficients[4], m4_d$coefficients[4],
           m1_r$coefficients[4], m2_r$coefficients[4], # reps coefs wokefiscal
           m3_r$coefficients[4], m4_r$coefficients[4])
se1 <- sqrt(diag(vcov(m1_d))); se2 <- sqrt(diag(vcov(m2_d))); se3 <- sqrt(diag(vcov(m3_d))); se4 <- sqrt(diag(vcov(m4_d))) # dems SEs 
se5 <- sqrt(diag(vcov(m1_r))); se6 <- sqrt(diag(vcov(m2_r))); se7 <- sqrt(diag(vcov(m3_r))); se8 <- sqrt(diag(vcov(m4_r))) # reps SEs 

hi1 <- c(coef1[1] +(1.96*se1[2]),   coef1[2]+(1.96*se2[2]),
         coef1[3] +(1.96*se3[2]),   coef1[4]+(1.96*se4[2]),
         coef1[5] +(1.96*se5[2]),   coef1[6]+(1.96*se6[2]),
         coef1[7] +(1.96*se7[2]),   coef1[8]+(1.96*se8[2]),
         coef1[9] +(1.96*se1[3]),  coef1[10]+(1.96*se2[3]),
         coef1[11]+(1.96*se3[3]),  coef1[12]+(1.96*se4[3]),
         coef1[13]+(1.96*se5[3]),  coef1[14]+(1.96*se6[3]),
         coef1[15]+(1.96*se7[3]),  coef1[16]+(1.96*se8[3]),
         coef1[17]+(1.96*se1[4]),  coef1[18]+(1.96*se2[4]),
         coef1[19]+(1.96*se3[4]),  coef1[20]+(1.96*se4[4]),
         coef1[21]+(1.96*se5[4]),  coef1[22]+(1.96*se6[4]),
         coef1[23]+(1.96*se7[4]),  coef1[24]+(1.96*se8[4]))

lo1 <- c(coef1[1] -(1.96*se1[2]),   coef1[2]-(1.96*se2[2]),
         coef1[3] -(1.96*se3[2]),   coef1[4]-(1.96*se4[2]),
         coef1[5] -(1.96*se5[2]),   coef1[6]-(1.96*se6[2]),
         coef1[7] -(1.96*se7[2]),   coef1[8]-(1.96*se8[2]),
         coef1[9] -(1.96*se1[3]),  coef1[10]-(1.96*se2[3]),
         coef1[11]-(1.96*se3[3]),  coef1[12]-(1.96*se4[3]),
         coef1[13]-(1.96*se5[3]),  coef1[14]-(1.96*se6[3]),
         coef1[15]-(1.96*se7[3]),  coef1[16]-(1.96*se8[3]),
         coef1[17]-(1.96*se1[4]),  coef1[18]-(1.96*se2[4]),
         coef1[19]-(1.96*se3[4]),  coef1[20]-(1.96*se4[4]),
         coef1[21]-(1.96*se5[4]),  coef1[22]-(1.96*se6[4]),
         coef1[23]-(1.96*se7[4]),  coef1[24]-(1.96*se8[4]))

party <- c('Democratic', 'Democratic', 'Democratic', 'Democratic',
           'Republican', 'Republican', 'Republican', 'Republican',
           'Democratic', 'Democratic', 'Democratic', 'Democratic',
           'Republican', 'Republican', 'Republican', 'Republican',
           'Democratic', 'Democratic', 'Democratic', 'Democratic',
           'Republican', 'Republican', 'Republican', 'Republican')
frame <- c('Fiscal', 'Fiscal', 'Fiscal', 'Fiscal',
           'Fiscal', 'Fiscal', 'Fiscal', 'Fiscal',
           'Woke', 'Woke', 'Woke', 'Woke',
           'Woke', 'Woke', 'Woke', 'Woke',
           'Woke Fiscal', 'Woke Fiscal', 'Woke Fiscal', 'Woke Fiscal',
           'Woke Fiscal', 'Woke Fiscal', 'Woke Fiscal', 'Woke Fiscal')
dv <- rep(c('CPT--Climate', 'CFS', 'CPT', 'Punish--Reward'), 6)

# Study I

coef1_c <- c(m1_d_c$coefficients[2], m2_d_c$coefficients[2], # dems coefs fiscal
           m3_d_c$coefficients[2], m4_d_c$coefficients[2],
           m1_r_c$coefficients[2], m2_r_c$coefficients[2], # reps coefs fiscal
           m3_r_c$coefficients[2], m4_r_c$coefficients[2],
           m1_d_c$coefficients[3], m2_d_c$coefficients[3], # dems coefs woke
           m3_d_c$coefficients[3], m4_d_c$coefficients[3],
           m1_r_c$coefficients[3], m2_r_c$coefficients[3], # reps coefs woke
           m3_r_c$coefficients[3], m4_r_c$coefficients[3],
           m1_d_c$coefficients[4], m2_d_c$coefficients[4], # dems coefs wokefiscal
           m3_d_c$coefficients[4], m4_d_c$coefficients[4],
           m1_r_c$coefficients[4], m2_r_c$coefficients[4], # reps coefs wokefiscal
           m3_r_c$coefficients[4], m4_r_c$coefficients[4])
se1_c <- sqrt(diag(vcov(m1_d_c))); se2_c <- sqrt(diag(vcov(m2_d_c))); se3_c <- sqrt(diag(vcov(m3_d_c))); se4_c <- sqrt(diag(vcov(m4_d_c))) # dems SEs 
se5_c <- sqrt(diag(vcov(m1_r_c))); se6_c <- sqrt(diag(vcov(m2_r_c))); se7_c <- sqrt(diag(vcov(m3_r_c))); se8_c <- sqrt(diag(vcov(m4_r_c))) # reps SEs 

hi1_c <- c(coef1_c[1] +(1.96*se1_c[2]),   coef1_c[2]+(1.96*se2_c[2]),
         coef1_c[3] +(1.96*se3_c[2]),   coef1_c[4]+(1.96*se4_c[2]),
         coef1_c[5] +(1.96*se5_c[2]),   coef1_c[6]+(1.96*se6_c[2]),
         coef1_c[7] +(1.96*se7_c[2]),   coef1_c[8]+(1.96*se8_c[2]),
         coef1_c[9] +(1.96*se1_c[3]),  coef1_c[10]+(1.96*se2_c[3]),
         coef1_c[11]+(1.96*se3_c[3]),  coef1_c[12]+(1.96*se4_c[3]),
         coef1_c[13]+(1.96*se5_c[3]),  coef1_c[14]+(1.96*se6_c[3]),
         coef1_c[15]+(1.96*se7_c[3]),  coef1_c[16]+(1.96*se8_c[3]),
         coef1_c[17]+(1.96*se1_c[4]),  coef1_c[18]+(1.96*se2_c[4]),
         coef1_c[19]+(1.96*se3_c[4]),  coef1_c[20]+(1.96*se4_c[4]),
         coef1_c[21]+(1.96*se5_c[4]),  coef1_c[22]+(1.96*se6_c[4]),
         coef1_c[23]+(1.96*se7_c[4]),  coef1_c[24]+(1.96*se8_c[4]))

lo1_c <- c(coef1_c[1] -(1.96*se1_c[2]),   coef1_c[2]-(1.96*se2_c[2]),
         coef1_c[3] -(1.96*se3_c[2]),   coef1_c[4]-(1.96*se4_c[2]),
         coef1_c[5] -(1.96*se5_c[2]),   coef1_c[6]-(1.96*se6_c[2]),
         coef1_c[7] -(1.96*se7_c[2]),   coef1_c[8]-(1.96*se8_c[2]),
         coef1_c[9] -(1.96*se1_c[3]),  coef1_c[10]-(1.96*se2_c[3]),
         coef1_c[11]-(1.96*se3_c[3]),  coef1_c[12]-(1.96*se4_c[3]),
         coef1_c[13]-(1.96*se5_c[3]),  coef1_c[14]-(1.96*se6_c[3]),
         coef1_c[15]-(1.96*se7_c[3]),  coef1_c[16]-(1.96*se8_c[3]),
         coef1_c[17]-(1.96*se1_c[4]),  coef1_c[18]-(1.96*se2_c[4]),
         coef1_c[19]-(1.96*se3_c[4]),  coef1_c[20]-(1.96*se4_c[4]),
         coef1_c[21]-(1.96*se5_c[4]),  coef1_c[22]-(1.96*se6_c[4]),
         coef1_c[23]-(1.96*se7_c[4]),  coef1_c[24]-(1.96*se8_c[4]))

party_c <- c('Democratic', 'Democratic', 'Democratic', 'Democratic',
           'Republican', 'Republican', 'Republican', 'Republican',
           'Democratic', 'Democratic', 'Democratic', 'Democratic',
           'Republican', 'Republican', 'Republican', 'Republican',
           'Democratic', 'Democratic', 'Democratic', 'Democratic',
           'Republican', 'Republican', 'Republican', 'Republican')
frame_c <- c('Fiscal', 'Fiscal', 'Fiscal', 'Fiscal',
           'Fiscal', 'Fiscal', 'Fiscal', 'Fiscal',
           'Woke', 'Woke', 'Woke', 'Woke',
           'Woke', 'Woke', 'Woke', 'Woke',
           'Woke Fiscal', 'Woke Fiscal', 'Woke Fiscal', 'Woke Fiscal',
           'Woke Fiscal', 'Woke Fiscal', 'Woke Fiscal', 'Woke Fiscal')
dv_c <- rep(c('CPT--Climate', 'CFS', 'CPT', 'Punish--Reward'), 6)


# Study II

coef2 <- c(m1_d2$coefficients[2], m2_d2$coefficients[2], # dems coefs fiscal
           m3_d2$coefficients[2], m4_d2$coefficients[2],
           m1_r2$coefficients[2], m2_r2$coefficients[2], # reps coefs fiscal
           m3_r2$coefficients[2], m4_r2$coefficients[2],
           m1_d2$coefficients[3], m2_d2$coefficients[3], # dems coefs woke
           m3_d2$coefficients[3], m4_d2$coefficients[3],
           m1_r2$coefficients[3], m2_r2$coefficients[3], # reps coefs woke
           m3_r2$coefficients[3], m4_r2$coefficients[3],
           m1_d2$coefficients[4], m2_d2$coefficients[4], # dems coefs wokefiscal
           m3_d2$coefficients[4], m4_d2$coefficients[4],
           m1_r2$coefficients[4], m2_r2$coefficients[4], # reps coefs wokefiscal
           m3_r2$coefficients[4], m4_r2$coefficients[4])
se1_2 <- sqrt(diag(vcov(m1_d2))); se2_2 <- sqrt(diag(vcov(m2_d2))); se3_2 <- sqrt(diag(vcov(m3_d2))); se4_2 <- sqrt(diag(vcov(m4_d2))) # dems SEs 
se5_2 <- sqrt(diag(vcov(m1_r2))); se6_2 <- sqrt(diag(vcov(m2_r2))); se7_2 <- sqrt(diag(vcov(m3_r2))); se8_2 <- sqrt(diag(vcov(m4_r2))) # reps SEs 

hi2 <- c(coef2[1] +(1.96*se1_2[2]),   coef2[2]+(1.96*se2_2[2]),
         coef2[3] +(1.96*se3_2[2]),   coef2[4]+(1.96*se4_2[2]),
         coef2[5] +(1.96*se5_2[2]),   coef2[6]+(1.96*se6_2[2]),
         coef2[7] +(1.96*se7_2[2]),   coef2[8]+(1.96*se8_2[2]),
         coef2[9] +(1.96*se1_2[3]),  coef2[10]+(1.96*se2_2[3]),
         coef2[11]+(1.96*se3_2[3]),  coef2[12]+(1.96*se4_2[3]),
         coef2[13]+(1.96*se5_2[3]),  coef2[14]+(1.96*se6_2[3]),
         coef2[15]+(1.96*se7_2[3]),  coef2[16]+(1.96*se8_2[3]),
         coef2[17]+(1.96*se1_2[4]),  coef2[18]+(1.96*se2_2[4]),
         coef2[19]+(1.96*se3_2[4]),  coef2[20]+(1.96*se4_2[4]),
         coef2[21]+(1.96*se5_2[4]),  coef2[22]+(1.96*se6_2[4]),
         coef2[23]+(1.96*se7_2[4]),  coef2[24]+(1.96*se8_2[4]))

lo2 <- c(coef2[1] -(1.96*se1_2[2]),   coef2[2]-(1.96*se2_2[2]),
         coef2[3] -(1.96*se3_2[2]),   coef2[4]-(1.96*se4_2[2]),
         coef2[5] -(1.96*se5_2[2]),   coef2[6]-(1.96*se6_2[2]),
         coef2[7] -(1.96*se7_2[2]),   coef2[8]-(1.96*se8_2[2]),
         coef2[9] -(1.96*se1_2[3]),  coef2[10]-(1.96*se2_2[3]),
         coef2[11]-(1.96*se3_2[3]),  coef2[12]-(1.96*se4_2[3]),
         coef2[13]-(1.96*se5_2[3]),  coef2[14]-(1.96*se6_2[3]),
         coef2[15]-(1.96*se7_2[3]),  coef2[16]-(1.96*se8_2[3]),
         coef2[17]-(1.96*se1_2[4]),  coef2[18]-(1.96*se2_2[4]),
         coef2[19]-(1.96*se3_2[4]),  coef2[20]-(1.96*se4_2[4]),
         coef2[21]-(1.96*se5_2[4]),  coef2[22]-(1.96*se6_2[4]),
         coef2[23]-(1.96*se7_2[4]),  coef2[24]-(1.96*se8_2[4]))

# party <- c('Democratic', 'Democratic', 'Democratic', 'Democratic',
#            'Republican', 'Republican', 'Republican', 'Republican',
#            'Democratic', 'Democratic', 'Democratic', 'Democratic',
#            'Republican', 'Republican', 'Republican', 'Republican',
#            'Democratic', 'Democratic', 'Democratic', 'Democratic',
#            'Republican', 'Republican', 'Republican', 'Republican')
# frame <- c('Fiscal', 'Fiscal', 'Fiscal', 'Fiscal',
#            'Fiscal', 'Fiscal', 'Fiscal', 'Fiscal',
#            'Woke', 'Woke', 'Woke', 'Woke',
#            'Woke', 'Woke', 'Woke', 'Woke',
#            'Woke Fiscal', 'Woke Fiscal', 'Woke Fiscal', 'Woke Fiscal',
#            'Woke Fiscal', 'Woke Fiscal', 'Woke Fiscal', 'Woke Fiscal')
# dv <- rep(c('CPT--Climate', 'CFS', 'CPT', 'Punish--Reward'), 6)

# Study 2 w/ controls

# Study II

coef2_c <- c(m1_d2_c$coefficients[2], m2_d2_c$coefficients[2], # dems coefs fiscal
           m3_d2_c$coefficients[2], m4_d2_c$coefficients[2],
           m1_r2_c$coefficients[2], m2_r2_c$coefficients[2], # reps coefs fiscal
           m3_r2_c$coefficients[2], m4_r2_c$coefficients[2],
           m1_d2_c$coefficients[3], m2_d2_c$coefficients[3], # dems coefs woke
           m3_d2_c$coefficients[3], m4_d2_c$coefficients[3],
           m1_r2_c$coefficients[3], m2_r2_c$coefficients[3], # reps coefs woke
           m3_r2_c$coefficients[3], m4_r2_c$coefficients[3],
           m1_d2_c$coefficients[4], m2_d2_c$coefficients[4], # dems coefs wokefiscal
           m3_d2_c$coefficients[4], m4_d2_c$coefficients[4],
           m1_r2_c$coefficients[4], m2_r2_c$coefficients[4], # reps coefs wokefiscal
           m3_r2_c$coefficients[4], m4_r2_c$coefficients[4])
se1_2_c <- sqrt(diag(vcov(m1_d2_c))); se2_2_c <- sqrt(diag(vcov(m2_d2_c))); se3_2_c <- sqrt(diag(vcov(m3_d2_c))); se4_2_c <- sqrt(diag(vcov(m4_d2_c))) # dems SEs 
se5_2_c <- sqrt(diag(vcov(m1_r2_c))); se6_2_c <- sqrt(diag(vcov(m2_r2_c))); se7_2_c <- sqrt(diag(vcov(m3_r2_c))); se8_2_c <- sqrt(diag(vcov(m4_r2_c))) # reps SEs 

hi2_c <- c(coef2_c[1] +(1.96*se1_2_c[2]),   coef2_c[2]+(1.96*se2_2_c[2]),
         coef2_c[3] +(1.96*se3_2_c[2]),   coef2_c[4]+(1.96*se4_2_c[2]),
         coef2_c[5] +(1.96*se5_2_c[2]),   coef2_c[6]+(1.96*se6_2_c[2]),
         coef2_c[7] +(1.96*se7_2_c[2]),   coef2_c[8]+(1.96*se8_2_c[2]),
         coef2_c[9] +(1.96*se1_2_c[3]),  coef2_c[10]+(1.96*se2_2_c[3]),
         coef2_c[11]+(1.96*se3_2_c[3]),  coef2_c[12]+(1.96*se4_2_c[3]),
         coef2_c[13]+(1.96*se5_2_c[3]),  coef2_c[14]+(1.96*se6_2_c[3]),
         coef2_c[15]+(1.96*se7_2_c[3]),  coef2_c[16]+(1.96*se8_2_c[3]),
         coef2_c[17]+(1.96*se1_2_c[4]),  coef2_c[18]+(1.96*se2_2_c[4]),
         coef2_c[19]+(1.96*se3_2_c[4]),  coef2_c[20]+(1.96*se4_2_c[4]),
         coef2_c[21]+(1.96*se5_2_c[4]),  coef2_c[22]+(1.96*se6_2_c[4]),
         coef2_c[23]+(1.96*se7_2_c[4]),  coef2_c[24]+(1.96*se8_2_c[4]))

lo2_c <- c(coef2_c[1] -(1.96*se1_2_c[2]),   coef2_c[2]-(1.96*se2_2_c[2]),
         coef2_c[3] -(1.96*se3_2_c[2]),   coef2_c[4]-(1.96*se4_2_c[2]),
         coef2_c[5] -(1.96*se5_2_c[2]),   coef2_c[6]-(1.96*se6_2_c[2]),
         coef2_c[7] -(1.96*se7_2_c[2]),   coef2_c[8]-(1.96*se8_2_c[2]),
         coef2_c[9] -(1.96*se1_2_c[3]),  coef2_c[10]-(1.96*se2_2_c[3]),
         coef2_c[11]-(1.96*se3_2_c[3]),  coef2_c[12]-(1.96*se4_2_c[3]),
         coef2_c[13]-(1.96*se5_2_c[3]),  coef2_c[14]-(1.96*se6_2_c[3]),
         coef2_c[15]-(1.96*se7_2_c[3]),  coef2_c[16]-(1.96*se8_2_c[3]),
         coef2_c[17]-(1.96*se1_2_c[4]),  coef2_c[18]-(1.96*se2_2_c[4]),
         coef2_c[19]-(1.96*se3_2_c[4]),  coef2_c[20]-(1.96*se4_2_c[4]),
         coef2_c[21]-(1.96*se5_2_c[4]),  coef2_c[22]-(1.96*se6_2_c[4]),
         coef2_c[23]-(1.96*se7_2_c[4]),  coef2_c[24]-(1.96*se8_2_c[4]))

study1 <- rep('Study I', 24)
study2 <- rep('Study II', 24)

df_1 <- as.data.frame(cbind(coef1, hi1, lo1, party, frame, dv, study1))
df_2 <- as.data.frame(cbind(coef2, hi2, lo2, party, frame, dv, study2))
names(df_1) <- c('coef', 'hi', 'lo', 'Party', 'frame', 'dv', 'study')
names(df_2) <- c('coef', 'hi', 'lo', 'Party', 'frame', 'dv', 'study')

df_1_c <- as.data.frame(cbind(coef1_c, hi1_c, lo1_c, party, frame, dv, study1))
df_2_c <- as.data.frame(cbind(coef2_c, hi2_c, lo2_c, party, frame, dv, study2))
names(df_1_c) <- c('coef', 'hi', 'lo', 'Party', 'frame', 'dv', 'study')
names(df_2_c) <- c('coef', 'hi', 'lo', 'Party', 'frame', 'dv', 'study')

df_1[, 1:3] <- lapply(df_1[, 1:3], as.numeric)
df_1[, 4:7] <- lapply(df_1[, 4:7], as.factor)
df_2[, 1:3] <- lapply(df_2[, 1:3], as.numeric)
df_2[, 4:7] <- lapply(df_2[, 4:7], as.factor)

df_1_c[, 1:3] <- lapply(df_1_c[, 1:3], as.numeric)
df_1_c[, 4:7] <- lapply(df_1_c[, 4:7], as.factor)
df_2_c[, 1:3] <- lapply(df_2_c[, 1:3], as.numeric)
df_2_c[, 4:7] <- lapply(df_2_c[, 4:7], as.factor)

df <- rbind(df_1, df_2)

p1 <- ggplot(df_1, aes(shape = dv, colour = Party)) +
  theme_bw() +
  geom_vline(xintercept = 0, linetype = 'dashed') +
  geom_pointrange(aes(x = coef, y = frame, xmin = lo,
                         xmax = hi),
                     lwd = 1, position = position_dodge(width = 0.5), fill = 'WHITE') +
  xlim(-0.2, 0.9) +
  facet_wrap(~Party) +
  scale_shape_manual(name = 'Dependent Variable',
                    labels = c('CFS', 'CPT', 'CPT-Climate', 'Punish-Reward'),
                    values = c(1, 17, 12, 19),
                    guide = guide_legend(reverse = T)) +
  scale_y_discrete(name = 'Treatment',
                   labels = c('Fiscal',
                              'Woke', 
                              'Woke Fiscal')) +
  scale_colour_manual(values = c(('darkblue'), "darkred")) +
  scale_x_continuous("Estimated Effect") +
  labs(title = 'Support for Corporate Political Engagement by Framing, Study I', caption = '\nNote: Estimated effects were calculated from two separate models (based on partisan identification). \nError bars represent 95% confidence intervals. \nFull results tables in Supplementary Information.') +
  theme(axis.text=element_text(size = 10)) + 
  theme(plot.title = element_text(hjust = 0.5)) + 
  theme(plot.caption = element_text(hjust = 0)) +
  theme(text=element_text(size = 10, family="serif"))

p2 <- ggplot(df_2, aes(shape = dv, colour = Party)) +
  theme_bw() +
  geom_vline(xintercept = 0, linetype = 'dashed') +
  geom_pointrange(aes(x = coef, y = frame, xmin = lo,
                         xmax = hi),
                     lwd = 1, position = position_dodge(width = 0.5), fill = 'WHITE') +
  xlim(-0.2, 0.9) +
  facet_wrap(~Party) +
  scale_shape_manual(name = 'Dependent Variable',
                    labels = c('CFS', 'CPT', 'CPT-Climate', 'Punish-Reward'),
                    values = c(1, 17, 12, 19),
                    guide = guide_legend(reverse = T)) +
  scale_y_discrete(name = 'Treatment',
                   labels = c('Fiscal',
                              'Woke', 
                              'Woke Fiscal')) +
  scale_colour_manual(values = c(('darkblue'), "darkred")) +
  scale_x_continuous("Estimated Effect") +
  labs(title = 'Support for Corporate Political Engagement by Framing, Study II', caption = '\nNote: Estimated effects were calculated from two separate models (based on partisan identification). \nResults estimated using ACS-based weights provided by YouGov. Error bars represent 95% confidence intervals. \nFull results tables in Supplementary Information.') +
  theme(axis.text=element_text(size = 10)) + 
  theme(plot.title = element_text(hjust = 0.5)) + 
  theme(plot.caption = element_text(hjust = 0)) +
  theme(text=element_text(size = 10, family="serif"))

p1; p2

# ggsave('cpt2_study1.png', p1, 'png', '~/Dropbox/Wayde/Research/Apps/Overleaf/CPT2_Woke/figures_tables_inputs',
#     width = 190, height = 127, units = 'mm')
# 
# ggsave('cpt2_study2.png', p2, 'png', '~/Dropbox/Wayde/Research/Apps/Overleaf/CPT2_Woke/figures_tables_inputs',
#     width = 190, height = 127, units = 'mm')

p1_c <- ggplot(df_1_c, aes(shape = dv, colour = Party)) +
  theme_bw() +
  geom_vline(xintercept = 0, linetype = 'dashed') +
  geom_pointrange(aes(x = coef, y = frame, xmin = lo,
                         xmax = hi),
                     lwd = 1, position = position_dodge(width = 0.5), fill = 'WHITE') +
  xlim(-0.2, 0.9) +
  facet_wrap(~Party) +
  scale_shape_manual(name = 'Dependent Variable',
                    labels = c('CFS', 'CPT', 'CPT-Climate', 'Punish-Reward'),
                    values = c(1, 17, 12, 19),
                    guide = guide_legend(reverse = T)) +
  scale_y_discrete(name = 'Treatment',
                   labels = c('Fiscal',
                              'Woke', 
                              'Woke Fiscal')) +
  scale_colour_manual(values = c(('darkblue'), "darkred")) +
  scale_x_continuous("Estimated Effect") +
  labs(title = 'Support for Corporate Political Engagement by Framing with Controls, Study I', caption = '\nNote: Estimated effects were calculated from two separate models (based on partisan identification). \nError bars represent 95% confidence intervals. Includes controls for age and highest level of education. \nFull results tables in Supplementary Information.') +
  theme(axis.text=element_text(size = 10)) + 
  theme(plot.title = element_text(hjust = 0.5)) + 
  theme(plot.caption = element_text(hjust = 0)) +
  theme(text=element_text(size = 10, family="serif"))

p2_c <- ggplot(df_2_c, aes(shape = dv, colour = Party)) +
  theme_bw() +
  geom_vline(xintercept = 0, linetype = 'dashed') +
  geom_pointrange(aes(x = coef, y = frame, xmin = lo,
                         xmax = hi),
                     lwd = 1, position = position_dodge(width = 0.5), fill = 'WHITE') +
  xlim(-0.2, 0.9) +
  facet_wrap(~Party) +
  scale_shape_manual(name = 'Dependent Variable',
                    labels = c('CFS', 'CPT', 'CPT-Climate', 'Punish-Reward'),
                    values = c(1, 17, 12, 19),
                    guide = guide_legend(reverse = T)) +
  scale_y_discrete(name = 'Treatment',
                   labels = c('Fiscal',
                              'Woke', 
                              'Woke Fiscal')) +
  scale_colour_manual(values = c(('darkblue'), "darkred")) +
  scale_x_continuous("Estimated Effect") +
  labs(title = 'Support for Corporate Political Engagement by Framing with Controls, Study II', caption = '\nNote: Estimated effects were calculated from two separate models (based on partisan identification). \nResults estimated using ACS-based weights provided by YouGov. Error bars represent 95% confidence intervals. Includes controls for age, highest level of education, and ideology. \nFull results tables in Supplementary Information.') +
  theme(axis.text=element_text(size = 10)) + 
  theme(plot.title = element_text(hjust = 0.5)) + 
  theme(plot.caption = element_text(hjust = 0)) +
  theme(text=element_text(size = 10, family="serif"))

p1_c; p2_c

# ggsave('cpt2_study1_c.png', p1_c, 'png', '~/Dropbox/Wayde/Research/Apps/Overleaf/CPT2_Woke/figures_tables_inputs',
#     width = 190, height = 127, units = 'mm')
# 
# ggsave('cpt2_study2_c.png', p2_c, 'png', '~/Dropbox/Wayde/Research/Apps/Overleaf/CPT2_Woke/figures_tables_inputs',
#     width = 190, height = 127, units = 'mm')

# Combined Figure

df[, 1:3] <- lapply(df[, 1:3], as.numeric)
df[, 4:7] <- lapply(df[, 4:7], as.factor)


plot1 <- ggplot(df, aes(shape = dv, colour = Party)) +
  theme_bw() +
  geom_vline(xintercept = 0, linetype = 'dashed') +
  geom_pointrange(aes(x = coef, y = frame, xmin = lo,
                         xmax = hi),
                     lwd = 1, position = position_dodge(width = 0.5), fill = 'WHITE') +
  xlim(-0.2, 0.9) +
  facet_wrap(~study) +
  scale_shape_manual(name = 'Dependent Variable',
                    labels = c('CFS', 'CPT', 'CPT-Climate', 'Punish-Reward'),
                    values = c(1, 17, 12, 19),
                    guide = guide_legend(reverse = T)) +
  scale_y_discrete(name = 'Treatment',
                   labels = c('Fiscal',
                              'Woke', 
                              'Woke Fiscal')) +
  scale_colour_manual(values = c(('darkblue'), "darkred")) +
  scale_x_continuous("Estimated Effect") +
  labs(title = 'Woke Framing Effects Support for Corporate Political Engagement', caption = '\nNote: Estimated effects were calculated from two separate models (based on partisan identification) for each study. \nResults from study II include ACS-based weights provided by YouGov. Error bars represent 95% confidence intervals. \nFull results tables in Supplementary Information.') +
  theme(axis.text=element_text(size = 10)) + 
  theme(plot.title = element_text(hjust = 0.5)) + 
  theme(plot.caption = element_text(hjust = 0)) +
  theme(text=element_text(size = 10, family="serif"))

plot1

# ggsave('cpt2_combo.png', plot1, 'png', '~/Dropbox/Wayde/Research/Apps/Overleaf/CPT2_Woke/figures_tables_inputs',
#     width = 190, height = 127, units = 'mm')

# ggsave('cpt2_combo.png', plot1, 'png', '/Users/wmarsh5/Dropbox/Wayde/Research/Apps/Overleaf/CPT2_Woke/figures_tables_inputs',
#     width = 190, height = 127, units = 'mm')

```

# Rainey TOST (90 CI)

```{r}

df_si1 <- data.frame(Estimate = c(clx(m1_r)[2, 1], # Fiscal
                                  clx(m1_d)[2, 1],
                                  clx(m2_r)[2, 1],
                                  clx(m2_d)[2, 1],
                                  clx(m3_r)[2, 1],
                                  clx(m3_d)[2, 1],
                                  clx(m4_r)[2, 1],
                                  clx(m4_d)[2, 1],
                     
                                  clx(m1_r)[3, 1], # Woke
                                  clx(m1_d)[3, 1],
                                  clx(m2_r)[3, 1],
                                  clx(m2_d)[3, 1],
                                  clx(m3_r)[3, 1],
                                  clx(m3_d)[3, 1],
                                  clx(m4_r)[3, 1],
                                  clx(m4_d)[3, 1],
                     
                                  clx(m1_r)[4, 1], # WokeFiscal
                                  clx(m1_d)[4, 1],
                                  clx(m2_r)[4, 1],
                                  clx(m2_d)[4, 1],
                                  clx(m3_r)[4, 1],
                                  clx(m3_d)[4, 1],
                                  clx(m4_r)[4, 1],
                                  clx(m4_d)[4, 1]),
                  lo = c(clx(m1_r)[2,1] - 1.64*clx(m1_r)[2, 2],
                         clx(m1_d)[2,1] - 1.64*clx(m1_d)[2, 2],
                         clx(m2_r)[2,1] - 1.64*clx(m2_r)[2, 2],
                         clx(m2_d)[2,1] - 1.64*clx(m2_d)[2, 2],
                         clx(m3_r)[2,1] - 1.64*clx(m3_r)[2, 2],
                         clx(m3_d)[2,1] - 1.64*clx(m3_d)[2, 2],
                         clx(m4_r)[2,1] - 1.64*clx(m4_r)[2, 2],
                         clx(m4_d)[2,1] - 1.64*clx(m4_d)[2, 2],
                         
                         clx(m1_r)[3,1] - 1.64*clx(m1_r)[3, 2],
                         clx(m1_d)[3,1] - 1.64*clx(m1_d)[3, 2],
                         clx(m2_r)[3,1] - 1.64*clx(m2_r)[3, 2],
                         clx(m2_d)[3,1] - 1.64*clx(m2_d)[3, 2],
                         clx(m3_r)[3,1] - 1.64*clx(m3_r)[3, 2],
                         clx(m3_d)[3,1] - 1.64*clx(m3_d)[3, 2],
                         clx(m4_r)[3,1] - 1.64*clx(m4_r)[3, 2],
                         clx(m4_d)[3,1] - 1.64*clx(m4_d)[3, 2],
                         
                         clx(m1_r)[4,1] - 1.64*clx(m1_r)[4, 2],
                         clx(m1_d)[4,1] - 1.64*clx(m1_d)[4, 2],
                         clx(m2_r)[4,1] - 1.64*clx(m2_r)[4, 2],
                         clx(m2_d)[4,1] - 1.64*clx(m2_d)[4, 2],
                         clx(m3_r)[4,1] - 1.64*clx(m3_r)[4, 2],
                         clx(m3_d)[4,1] - 1.64*clx(m3_d)[4, 2],
                         clx(m4_r)[4,1] - 1.64*clx(m4_r)[4, 2],
                         clx(m4_d)[4,1] - 1.64*clx(m4_d)[4, 2]),
                  hi = c(clx(m1_r)[2,1] + 1.64*clx(m1_r)[2, 2],
                         clx(m1_d)[2,1] + 1.64*clx(m1_d)[2, 2],
                         clx(m2_r)[2,1] + 1.64*clx(m2_r)[2, 2],
                         clx(m2_d)[2,1] + 1.64*clx(m2_d)[2, 2],
                         clx(m3_r)[2,1] + 1.64*clx(m3_r)[2, 2],
                         clx(m3_d)[2,1] + 1.64*clx(m3_d)[2, 2],
                         clx(m4_r)[2,1] + 1.64*clx(m4_r)[2, 2],
                         clx(m4_d)[2,1] + 1.64*clx(m4_d)[2, 2],
                         
                         clx(m1_r)[3,1] + 1.64*clx(m1_r)[3, 2],
                         clx(m1_d)[3,1] + 1.64*clx(m1_d)[3, 2],
                         clx(m2_r)[3,1] + 1.64*clx(m2_r)[3, 2],
                         clx(m2_d)[3,1] + 1.64*clx(m2_d)[3, 2],
                         clx(m3_r)[3,1] + 1.64*clx(m3_r)[3, 2],
                         clx(m3_d)[3,1] + 1.64*clx(m3_d)[3, 2],
                         clx(m4_r)[3,1] + 1.64*clx(m4_r)[3, 2],
                         clx(m4_d)[3,1] + 1.64*clx(m4_d)[3, 2],
                         
                         clx(m1_r)[4,1] + 1.64*clx(m1_r)[4, 2],
                         clx(m1_d)[4,1] + 1.64*clx(m1_d)[4, 2],
                         clx(m2_r)[4,1] + 1.64*clx(m2_r)[4, 2],
                         clx(m2_d)[4,1] + 1.64*clx(m2_d)[4, 2],
                         clx(m3_r)[4,1] + 1.64*clx(m3_r)[4, 2],
                         clx(m3_d)[4,1] + 1.64*clx(m3_d)[4, 2],
                         clx(m4_r)[4,1] + 1.64*clx(m4_r)[4, 2],
                         clx(m4_d)[4,1] + 1.64*clx(m4_d)[4, 2]),
                  treatment = c(rep('Fiscal', 8), rep('Woke', 8), rep('WokeFiscal', 8)),
                  DV = rep(c('CPT-Climate', 'CPT-Climate', 'CFS', 'CFS', 'CPT', 'CPT', 'Punish-Reward', 'Punish-Reward'), 3),
                  Party = rep(c('Republican', 'Democrat'), 12))

df_si2 <- data.frame(Estimate = c(clx(m1_r2)[2, 1], # Fiscal
                                  clx(m1_d2)[2, 1],
                                  clx(m2_r2)[2, 1],
                                  clx(m2_d2)[2, 1],
                                  clx(m3_r2)[2, 1],
                                  clx(m3_d2)[2, 1],
                                  clx(m4_r2)[2, 1],
                                  clx(m4_d2)[2, 1],
                     
                                  clx(m1_r2)[3, 1], # Woke
                                  clx(m1_d2)[3, 1],
                                  clx(m2_r2)[3, 1],
                                  clx(m2_d2)[3, 1],
                                  clx(m3_r2)[3, 1],
                                  clx(m3_d2)[3, 1],
                                  clx(m4_r2)[3, 1],
                                  clx(m4_d2)[3, 1],
                     
                                  clx(m1_r2)[4, 1], # WokeFiscal
                                  clx(m1_d2)[4, 1],
                                  clx(m2_r2)[4, 1],
                                  clx(m2_d2)[4, 1],
                                  clx(m3_r2)[4, 1],
                                  clx(m3_d2)[4, 1],
                                  clx(m4_r2)[4, 1],
                                  clx(m4_d2)[4, 1]),
                  lo = c(clx(m1_r2)[2,1] - 1.64*clx(m1_r2)[2, 2],
                         clx(m1_d2)[2,1] - 1.64*clx(m1_d2)[2, 2],
                         clx(m2_r2)[2,1] - 1.64*clx(m2_r2)[2, 2],
                         clx(m2_d2)[2,1] - 1.64*clx(m2_d2)[2, 2],
                         clx(m3_r2)[2,1] - 1.64*clx(m3_r2)[2, 2],
                         clx(m3_d2)[2,1] - 1.64*clx(m3_d2)[2, 2],
                         clx(m4_r2)[2,1] - 1.64*clx(m4_r2)[2, 2],
                         clx(m4_d2)[2,1] - 1.64*clx(m4_d2)[2, 2],
                         
                         clx(m1_r2)[3,1] - 1.64*clx(m1_r2)[3, 2],
                         clx(m1_d2)[3,1] - 1.64*clx(m1_d2)[3, 2],
                         clx(m2_r2)[3,1] - 1.64*clx(m2_r2)[3, 2],
                         clx(m2_d2)[3,1] - 1.64*clx(m2_d2)[3, 2],
                         clx(m3_r2)[3,1] - 1.64*clx(m3_r2)[3, 2],
                         clx(m3_d2)[3,1] - 1.64*clx(m3_d2)[3, 2],
                         clx(m4_r2)[3,1] - 1.64*clx(m4_r2)[3, 2],
                         clx(m4_d2)[3,1] - 1.64*clx(m4_d2)[3, 2],
                         
                         clx(m1_r2)[4,1] - 1.64*clx(m1_r2)[4, 2],
                         clx(m1_d2)[4,1] - 1.64*clx(m1_d2)[4, 2],
                         clx(m2_r2)[4,1] - 1.64*clx(m2_r2)[4, 2],
                         clx(m2_d2)[4,1] - 1.64*clx(m2_d2)[4, 2],
                         clx(m3_r2)[4,1] - 1.64*clx(m3_r2)[4, 2],
                         clx(m3_d2)[4,1] - 1.64*clx(m3_d2)[4, 2],
                         clx(m4_r2)[4,1] - 1.64*clx(m4_r2)[4, 2],
                         clx(m4_d2)[4,1] - 1.64*clx(m4_d2)[4, 2]),
                  hi = c(clx(m1_r2)[2,1] + 1.64*clx(m1_r2)[2, 2],
                         clx(m1_d2)[2,1] + 1.64*clx(m1_d2)[2, 2],
                         clx(m2_r2)[2,1] + 1.64*clx(m2_r2)[2, 2],
                         clx(m2_d2)[2,1] + 1.64*clx(m2_d2)[2, 2],
                         clx(m3_r2)[2,1] + 1.64*clx(m3_r2)[2, 2],
                         clx(m3_d2)[2,1] + 1.64*clx(m3_d2)[2, 2],
                         clx(m4_r2)[2,1] + 1.64*clx(m4_r2)[2, 2],
                         clx(m4_d2)[2,1] + 1.64*clx(m4_d2)[2, 2],
                         
                         clx(m1_r2)[3,1] + 1.64*clx(m1_r2)[3, 2],
                         clx(m1_d2)[3,1] + 1.64*clx(m1_d2)[3, 2],
                         clx(m2_r2)[3,1] + 1.64*clx(m2_r2)[3, 2],
                         clx(m2_d2)[3,1] + 1.64*clx(m2_d2)[3, 2],
                         clx(m3_r2)[3,1] + 1.64*clx(m3_r2)[3, 2],
                         clx(m3_d2)[3,1] + 1.64*clx(m3_d2)[3, 2],
                         clx(m4_r2)[3,1] + 1.64*clx(m4_r2)[3, 2],
                         clx(m4_d2)[3,1] + 1.64*clx(m4_d2)[3, 2],
                         
                         clx(m1_r2)[4,1] + 1.64*clx(m1_r2)[4, 2],
                         clx(m1_d2)[4,1] + 1.64*clx(m1_d2)[4, 2],
                         clx(m2_r2)[4,1] + 1.64*clx(m2_r2)[4, 2],
                         clx(m2_d2)[4,1] + 1.64*clx(m2_d2)[4, 2],
                         clx(m3_r2)[4,1] + 1.64*clx(m3_r2)[4, 2],
                         clx(m3_d2)[4,1] + 1.64*clx(m3_d2)[4, 2],
                         clx(m4_r2)[4,1] + 1.64*clx(m4_r2)[4, 2],
                         clx(m4_d2)[4,1] + 1.64*clx(m4_d2)[4, 2]),
                  treatment = c(rep('Fiscal', 8), rep('Woke', 8), rep('WokeFiscal', 8)),
                  DV = rep(c('CPT-Climate', 'CPT-Climate', 'CFS', 'CFS', 'CPT', 'CPT', 'Punish-Reward', 'Punish-Reward'), 3),
                  Party = rep(c('Republican', 'Democrat'), 12))

df_si <- rbind(df_si1, df_si2)
df_si$Study <- c(rep('Study I', 24), rep('Study II', 24))

df_si_r <- subset(df_si, Party == 'Republican')
df_si_d <- subset(df_si, Party == 'Democrat')

m <- 0.25

fig_si_r <- ggplot(df_si_r, aes(shape = dv)) + 
  theme_bw() +
  coord_flip() +
  facet_wrap(~Study) +
  scale_shape_manual(name = 'Dependent Variable',
                    labels = c('CFS', 'CPT', 'CPT-Climate', 'Punish-Reward'),
                    values = c(1, 17, 12, 19),
                    guide = guide_legend(reverse = T)) +
  geom_vline(xintercept = 0, linetype = 'solid') +
  geom_vline(xintercept = c(-m, m), linetype = "dashed") + 
  geom_pointrange(aes(x = Estimate, y = treatment, xmin = lo,
                         xmax = hi),
                     lwd = 1, position = position_dodge(width = 0.5), fill = 'WHITE') +
  xlim(-1, 1) +
  xlab('Treatment Effect') + 
  ylab('') +
  labs(title = 'Rainey (2014) Test of Negligible Effects, Republican Respondents') +
  theme(axis.text=element_text(size = 10)) + 
  theme(plot.title = element_text(hjust = 0.5)) + 
  theme(plot.caption = element_text(hjust = 0)) +
  theme(text=element_text(size = 10, family="serif"))

fig_si_d <- ggplot(df_si_d, aes(shape = dv)) + 
  theme_bw() +
  coord_flip() +
  facet_wrap(~Study) +
  scale_shape_manual(name = 'Dependent Variable',
                    labels = c('CFS', 'CPT', 'CPT-Climate', 'Punish-Reward'),
                    values = c(1, 17, 12, 19),
                    guide = guide_legend(reverse = T)) + 
  geom_vline(xintercept = 0, linetype = 'solid') +
  geom_vline(xintercept = c(-m, m), linetype = "dashed") +
  geom_pointrange(aes(x = Estimate, y = treatment, xmin = lo,
                         xmax = hi),
                     lwd = 1, position = position_dodge(width = 0.5), fill = 'WHITE') +
  xlim(-1, 1) +
  xlab('Treatment Effect') + 
  ylab('') +
  labs(title = 'Rainey (2014) Test of Negligible Effects, Democratic Respondents') +
  theme(axis.text=element_text(size = 10)) + 
  theme(plot.title = element_text(hjust = 0.5)) + 
  theme(plot.caption = element_text(hjust = 0)) +
  theme(text=element_text(size = 10, family="serif"))


ggsave('rainey_d.png', fig_si_d, 'png', '~/Dropbox/Wayde/Research/Apps/Overleaf/CPT2_Woke/figures_tables_inputs',
    width = 190, height = 127, units = 'mm')

ggsave('rainey_r.png', fig_si_r, 'png', '~/Dropbox/Wayde/Research/Apps/Overleaf/CPT2_Woke/figures_tables_inputs',
    width = 190, height = 127, units = 'mm')

```

## Bonferroni-Corrected Pair-Wise T-Tests

```{r Bonferroni}

df_r <- subset(dat1, pid2 == 'Republican')
df_d <- subset(dat1, pid2 == 'Democrat')

df_r2 <- subset(dat2, pid2 == 'Republican')
df_d2 <- subset(dat2, pid2 == 'Democrat')

## Study I
bon1 <- pairwise.t.test(dat1$dv1, dat1$Frame, p.adjust.method="bonferroni")
bon2 <- pairwise.t.test(dat1$dv2, dat1$Frame, p.adjust.method="bonferroni")
bon3 <- pairwise.t.test(dat1$dv3, dat1$Frame, p.adjust.method="bonferroni")
bon4 <- pairwise.t.test(dat1$dv4, dat1$Frame, p.adjust.method="bonferroni")

bon1r <- pairwise.t.test(df_r$dv1, df_r$Frame, p.adjust.method="bonferroni")
bon2r <- pairwise.t.test(df_r$dv2, df_r$Frame, p.adjust.method="bonferroni")
bon3r <- pairwise.t.test(df_r$dv3, df_r$Frame, p.adjust.method="bonferroni")
bon4r <- pairwise.t.test(df_r$dv4, df_r$Frame, p.adjust.method="bonferroni")

bon1d <- pairwise.t.test(df_d$dv1, df_d$Frame, p.adjust.method="bonferroni")
bon2d <- pairwise.t.test(df_d$dv2, df_d$Frame, p.adjust.method="bonferroni")
bon3d <- pairwise.t.test(df_d$dv3, df_d$Frame, p.adjust.method="bonferroni")
bon4d <- pairwise.t.test(df_d$dv4, df_d$Frame, p.adjust.method="bonferroni")

## Study II
bon12 <- pairwise.t.test(dat2$dv1, dat2$Frame, p.adjust.method="bonferroni")
bon22 <- pairwise.t.test(dat2$dv2, dat2$Frame, p.adjust.method="bonferroni")
bon32 <- pairwise.t.test(dat2$dv3, dat2$Frame, p.adjust.method="bonferroni")
bon42 <- pairwise.t.test(dat2$dv4, dat2$Frame, p.adjust.method="bonferroni")

bon1r2 <- pairwise.t.test(df_r2$dv1, df_r2$Frame, p.adjust.method="bonferroni")
bon2r2 <- pairwise.t.test(df_r2$dv2, df_r2$Frame, p.adjust.method="bonferroni")
bon3r2 <- pairwise.t.test(df_r2$dv3, df_r2$Frame, p.adjust.method="bonferroni")
bon4r2 <- pairwise.t.test(df_r2$dv4, df_r2$Frame, p.adjust.method="bonferroni")

bon1d2 <- pairwise.t.test(df_d2$dv1, df_d2$Frame, p.adjust.method="bonferroni")
bon2d2 <- pairwise.t.test(df_d2$dv2, df_d2$Frame, p.adjust.method="bonferroni")
bon3d2 <- pairwise.t.test(df_d2$dv3, df_d2$Frame, p.adjust.method="bonferroni")
bon4d2 <- pairwise.t.test(df_d2$dv4, df_d2$Frame, p.adjust.method="bonferroni")

```

# Analysis of Exploratory DV

```{r Exploratory Analysis}

dat2$cpt_dv5 <- relevel(as.factor(dat1$cpt_dv5), ref = 'State legislatures should do nothing')
dat1_d <- subset(dat1, pid2 = 'Democrat')
dat1_r <- subset(dat1, pid2 = 'Republican')

m5_d <- multinom(cpt_dv5 ~ cpt_exp, data = dat1_d)
summary(m5_d)
m5_r <- multinom(cpt_dv5 ~ cpt_exp, data = dat1_r)
summary(m5_r)

p5_1 <- ggeffect(m5) %>%
  plot()



# m5_d <- model.matrix(cpt_dv5 ~ cpt_exp, data = dat1_d)
# m5_r <- model.matrix(cpt_dv5 ~ cpt_exp, data = dat1_r)
# 
# m5_d_coef <- m5_d%*%t(coef(m5_test_d))
# m5_r_coef <- m5_r%*%t(coef(m5_test_r))
# 
# m5_d_coef_exp <- exp(m5_d_coef)
# m5_r_coef_exp <- exp(m5_r_coef)
# pred.odd_d <- data.frame(Control = 1, m5_d_coef_exp)
# pred.odd_r <- data.frame(Control = 1, m5_r_coef_exp)
# 
# probmulti_d <- pred.odd_d/apply(pred.odd_d, 1, sum)
# probmulti_r <- pred.odd_r/apply(pred.odd_r, 1, sum)

```

## DIMs (Tukey Multiple pairwise-comparisons)

```{r}

aov_1r <- aov(cpt_dv1 ~ as.factor(cpt_exp), data = dat1_r)
dim_1r <- TukeyHSD(aov_1r)
aov_2r <- aov(cpt_dv2 ~ as.factor(cpt_exp), data = dat1_r)
dim_2r <- TukeyHSD(aov_2r)
aov_3r <- aov(cpt_dv3 ~ as.factor(cpt_exp), data = dat1_r)
dim_3r <- TukeyHSD(aov_3r)
aov_4r <- aov(cpt_dv4 ~ as.factor(cpt_exp), data = dat1_r)
dim_4r <- TukeyHSD(aov_4r)

aov_1d <- aov(cpt_dv1 ~ as.factor(cpt_exp), data = dat1_d)
dim_1d <- TukeyHSD(aov_1d)
aov_2d <- aov(cpt_dv2 ~ as.factor(cpt_exp), data = dat1_d)
dim_2d <- TukeyHSD(aov_2d)
aov_3d <- aov(cpt_dv3 ~ as.factor(cpt_exp), data = dat1_d)
dim_3d <- TukeyHSD(aov_3d)
aov_4d <- aov(cpt_dv4 ~ as.factor(cpt_exp), data = dat1_d)
dim_4d <- TukeyHSD(aov_4d)

dim_1r <- tidy(dim_1r)
dim_2r <- tidy(dim_2r)
dim_1d <- tidy(dim_1d)
dim_2d <- tidy(dim_2d)
dim_3r <- tidy(dim_3r)
dim_4r <- tidy(dim_4r)
dim_3d <- tidy(dim_3d)
dim_4d <- tidy(dim_4d)

library(broom); library(xtable)
print.xtable(xtable(dim_1r[c(2, 4, 7)]), include.rownames = F)
print.xtable(xtable(dim_1d[c(2, 4, 7)]), include.rownames = F)
print.xtable(xtable(dim_2r[c(2, 4, 7)]), include.rownames = F)
print.xtable(xtable(dim_2d[c(2, 4, 7)]), include.rownames = F)
print.xtable(xtable(dim_3r[c(2, 4, 7)]), include.rownames = F)
print.xtable(xtable(dim_3d[c(2, 4, 7)]), include.rownames = F)
print.xtable(xtable(dim_4r[c(2, 4, 7)]), include.rownames = F)
print.xtable(xtable(dim_4d[c(2, 4, 7)]), include.rownames = F)

## DV 5 figure

df_5 <- as.data.frame(rbind(dim_5d[1:3, ], dim_5r[1:3, ]))

df_5$pid <- c(rep('Democrat', 3), rep('Republican', 3))

plot5 <- ggplot(df_5, aes()) +
  theme_bw() +
  geom_vline(xintercept = 0, linetype = 'dashed') +
  geom_pointrange(aes(x = estimate, y = contrast, xmin = conf.low,
                         xmax = conf.high),
                     lwd = 1, position = position_dodge(width = 0.5), fill = 'WHITE') +
  xlim(-0.7, 0.7) +
  facet_wrap(~pid) +
  labs(x = 'Estimated Effect',
       y = 'Treatment Comparison')
  scale_y_discrete(name = 'Treatment',
                   labels = c('Fiscal',
                              'Woke', 
                              'Woke Fiscal')) +
  scale_x_continuous("Estimated Effect") +
  labs(title = 'Woke Framing Effects Support for \nCorporate Political Engagement', caption = '\nNote: Estimated effects were calculated from two separate models (based on partisan identification). \nError bars represent 95/% confidence intervals. Full results table in Supplementary Information.') +
  theme(axis.text=element_text(size = 10)) + 
  theme(plot.title = element_text(hjust = 0.5)) + 
  theme(plot.caption = element_text(hjust = 0)) +
  theme(text=element_text(size = 10, family="serif"))

plot5

```

# Summary Statistics and Diagnostics

## Survey Stats

```{r}

dat1_sums <- dat1[, c('gender', 'hispanic', 'race', 'age', 'educ', 'partyid')]

dat1_sums[1, 'age'] <- 30
dat1_sums[675, 'age'] <- 37

dat1_sums$female <- ifelse(dat1_sums$gender == 'Female', 1, 0)
dat1_sums$hispanic <- ifelse(dat1_sums$hispanic == 'Yes', 1, 0)
dat1_sums$college <- ifelse(dat1_sums$educ == "Bachelor's Degree" | 
                            dat1_sums$educ == "Postgraduate degree" , 1, 0)
dat1_sums$black <- ifelse(dat1_sums$race == 'Black', 1, 0)

dat1_sums$gender <- NULL
dat1_sums$race <- NULL
dat1_sums$educ <- NULL

round(colMeans(dat1_sums, na.rm = T), 2) 

dat2_sums <- dat2[, c('sex', 'hispanic', 'race', 'age', 'edu', 'pid7')]

dat2_sums$female <- ifelse(dat2_sums$sex > 0, 1, 0)
dat2_sums$hispanic <- ifelse(dat2_sums$hispanic == 1, 1, 0)
dat2_sums$college <- ifelse(dat2_sums$edu > 0.5, 1, 0)
dat2_sums$black <- ifelse(dat2_sums$race == 2, 1, 0)

dat2_sums$sex <- NULL
dat2_sums$race <- NULL
dat2_sums$edu <- NULL

round(colMeans(dat2_sums, na.rm = T), 2)

# Descriptives of DVs

hist.df <- as.data.frame(rbind(dat1[, c('cpt_dv1', 'cpt_dv2', 'cpt_dv3', 'cpt_dv4')]))
hist.df$study <- c(rep('Study I', 1521))
names(hist.df) <- c('CPT_Climate', 'CPS', 'CPT', 'Punish_Reward')


dv_hist1 <- ggplot(hist.df, aes(x = CPT_Climate)) + 
  geom_histogram(aes(y = ..density..), binwidth = 0.25) +
  scale_y_continuous(labels = scales::percent_format()) +
  # facet_wrap(~Study) +
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(legend.position="bottom") +
  theme(text=element_text(family="serif"))

dv_hist2 <- ggplot(hist.df, aes(x = CPS)) + 
  geom_histogram(aes(y = ..density..), binwidth = 0.25) +
  scale_y_continuous(labels = scales::percent_format()) +
  # facet_wrap(~Study) +
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(legend.position="bottom") +
  theme(text=element_text(family="serif"))

dv_hist3 <- ggplot(hist.df, aes(x = CPT)) + 
  geom_histogram(aes(y = ..density..), binwidth = 0.01) +
  scale_y_continuous(labels = scales::percent_format()) +
  # facet_wrap(~Study) +
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(legend.position="bottom") +
  theme(text=element_text(family="serif"))

dv_hist4 <- ggplot(hist.df, aes(x = Punish_Reward)) + 
  geom_histogram(aes(y = ..density..), binwidth = 0.5) +
  scale_y_continuous(labels = scales::percent_format()) +
  # facet_wrap(~Study) +
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(legend.position="bottom") +
  theme(text=element_text(family="serif"))

# ggsave('dv_hist1.png', dv_hist1, 'png', '~/Dropbox/Wayde/Research/Apps/Overleaf/CPT2_Woke/figures_tables_inputs',
#     width = 190, height = 127, units = 'mm')
# ggsave('dv_hist2.png', dv_hist2, 'png', '~/Dropbox/Wayde/Research/Apps/Overleaf/CPT2_Woke/figures_tables_inputs',
#     width = 190, height = 127, units = 'mm')
# ggsave('dv_hist3.png', dv_hist3, 'png', '~/Dropbox/Wayde/Research/Apps/Overleaf/CPT2_Woke/figures_tables_inputs',
#     width = 190, height = 127, units = 'mm')
# ggsave('dv_hist4.png', dv_hist4, 'png', '~/Dropbox/Wayde/Research/Apps/Overleaf/CPT2_Woke/figures_tables_inputs',
#     width = 190, height = 127, units = 'mm')

hist.df2 <- as.data.frame(rbind(dat2[, c('cpt_dv1', 'cpt_dv2', 'cpt_dv3', 'cpt_dv4')]))
hist.df2$study <- c(rep('Study II', 1000))
names(hist.df2) <- c('CPT_Climate', 'CPS', 'CPT', 'Punish_Reward')


dv_hist1_2 <- ggplot(hist.df2, aes(x = CPT_Climate)) + 
  geom_histogram(aes(y = ..density..), binwidth = 0.25) +
  scale_y_continuous(labels = scales::percent_format()) +
  # facet_wrap(~Study) +
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(legend.position="bottom") +
  theme(text=element_text(family="serif"))

dv_hist2_2 <- ggplot(hist.df2, aes(x = CPS)) + 
  geom_histogram(aes(y = ..density..), binwidth = 0.25) +
  scale_y_continuous(labels = scales::percent_format()) +
  # facet_wrap(~Study) +
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(legend.position="bottom") +
  theme(text=element_text(family="serif"))

dv_hist3_2 <- ggplot(hist.df2, aes(x = CPT)) + 
  geom_histogram(aes(y = ..density..), binwidth = 0.01) +
  scale_y_continuous(labels = scales::percent_format()) +
  # facet_wrap(~Study) +
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(legend.position="bottom") +
  theme(text=element_text(family="serif"))

dv_hist4_2 <- ggplot(hist.df2, aes(x = Punish_Reward)) + 
  geom_histogram(aes(y = ..density..), binwidth = 0.5) +
  scale_y_continuous(labels = scales::percent_format()) +
  # facet_wrap(~Study) +
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(legend.position="bottom") +
  theme(text=element_text(family="serif"))

ggsave('dv_hist1_2.png', dv_hist1_2, 'png', '~/Dropbox/Wayde/Research/Apps/Overleaf/CPT2_Woke/figures_tables_inputs',
    width = 190, height = 127, units = 'mm')
ggsave('dv_hist2_2.png', dv_hist2_2, 'png', '~/Dropbox/Wayde/Research/Apps/Overleaf/CPT2_Woke/figures_tables_inputs',
    width = 190, height = 127, units = 'mm')
ggsave('dv_hist3_2.png', dv_hist3_2, 'png', '~/Dropbox/Wayde/Research/Apps/Overleaf/CPT2_Woke/figures_tables_inputs',
    width = 190, height = 127, units = 'mm')
ggsave('dv_hist4_2.png', dv_hist4_2, 'png', '~/Dropbox/Wayde/Research/Apps/Overleaf/CPT2_Woke/figures_tables_inputs',
    width = 190, height = 127, units = 'mm')

```

## Post-Hoc Power Analysis

```{r}

power.f.ancova(eta.squared = 0.01,
               factor.levels = 4,
               power = 0.80,
               alpha = 0.05)

```

## Balance Analysis

```{r}

dat1_sums$tr1 <- dat1$cpt_exp
df_complete <- na.omit(dat1_sums)

model1_age <- aov(age ~ as.factor(tr1), data = df_complete)
model1_col <- aov(college ~ as.factor(tr1), data = df_complete)
model1_sex <- aov(female ~ as.factor(tr1), data = df_complete)
model1_blk <- aov(black ~ as.factor(tr1), data = df_complete)
model1_hsp <- aov(hispanic ~ as.factor(tr1), data = df_complete)
model1_pid <- aov(partyid ~ as.factor(tr1), data = df_complete)

gg_tukey<-function(Tukey){
  A<-require("tidyverse")
  if(A==TRUE){
    library(tidyverse)
  } else {
    install.packages("tidyverse")
    library(tidyverse)
  }
  B<-as.data.frame(Tukey[1])
  colnames(B)[2:3]<-c("min",
                      "max")
  C<-data.frame(id=row.names(B),
                min=B$min,
                max=B$max)
  D<-C%>%
    ggplot(aes(id))+
    geom_errorbar(aes(ymin=min,
                     ymax=max),
                 width = 0.2)+
    geom_hline(yintercept=0,
               color="red")+
    labs(x=NULL)+
    coord_flip()+
    theme(text=element_text(family="serif"),
            title=element_text(color="black",size=15),
            axis.text = element_text(color="black",size=10),
            axis.title = element_text(color="black",size=10),
            panel.grid=element_line(color="grey75"),
            axis.line=element_blank(),
            plot.background=element_rect(fill="white",color="white"),
            panel.background=element_rect(fill="white"),
            panel.border = element_rect(colour = "black", fill = NA,size=0.59),
            legend.key= element_rect(color="white",fill="white")
          )
  return(D)
}

tk1_1 <- gg_tukey(TukeyHSD(model1_age))
tk1_2 <- gg_tukey(TukeyHSD(model1_col))
tk1_3 <- gg_tukey(TukeyHSD(model1_sex))
tk1_4 <- gg_tukey(TukeyHSD(model1_blk))
tk1_5 <- gg_tukey(TukeyHSD(model1_hsp))
tk1_6 <- gg_tukey(TukeyHSD(model1_pid))

# tk1_1 <- plot(TukeyHSD(model1_age, conf.level = 0.95), las = 2)
# tk1_2 <- plot(TukeyHSD(model1_col, conf.level = 0.95), las = 2)
# tk1_3 <- plot(TukeyHSD(model1_sex, conf.level = 0.95), las = 2)
# tk1_4 <- plot(TukeyHSD(model1_blk, conf.level = 0.95), las = 2)
# tk1_5 <- plot(TukeyHSD(model1_hsp, conf.level = 0.95), las = 2)
# tk1_6 <- plot(TukeyHSD(model1_pid, conf.level = 0.95), las = 2)

tk1 <- ggarrange(tk1_1, tk1_2, tk1_3, tk1_4, tk1_5, tk1_6,
          labels = c('Age', 'College', 'Sex', 'Black', 'Hispanic', 'Party ID'),
          ncol = 2, nrow = 3)

```

# Causal Random Forests (from Green and White Elements)

```{r CRF Packages}

library(grf)
library(data.table)
library(mltools)
library(tidyverse)

```


```{r CRF Setup}

dat1_crf <- dat1[, c('gender', 'hispanic', 'race', 'age', 'educ', 'partyid', 'cpt_exp', 'cpt_dv1', 'cpt_dv2', 'cpt_dv3', 'cpt_dv4')]

dat1_crf[1, 'age'] <- 30
dat1_crf[675, 'age'] <- 37

dat1_crf$female <- ifelse(dat1_crf$gender == 'Female', 1, 0)
dat1_crf$hispanic <- ifelse(dat1_crf$hispanic == 'Yes', 1, 0)
dat1_crf$college <- ifelse(dat1_crf$educ == "Bachelor's Degree" | 
                            dat1_crf$educ == "Postgraduate degree" , 1, 0)
dat1_crf$black <- ifelse(dat1_crf$race == 'Black', 1, 0)
dat1_crf <- ifelse(dat1_crf$cpt_exp == 'Control')

dat1_crf$gender <- NULL
dat1_crf$race <- NULL
dat1_crf$educ <- NULL

data <- na.omit(dat1_crf)

# Define outcome and Treatments as binary vectors (should already be done)
# OUTCOME = openletter (or maybe a new variable checking what they wrote, tr_letter); ft_imm; immgood; increaseimm; refugeesupport; donation (or maybe checking to see amount  > 0 and call tr_donation)
# TREATMENT = FH vs. PT vs. Control, so probably two treatments tr_fh (1 if FH, 0 if control) and tr_pt (1 if PT, 0 if control)

Y <- as.numeric(dat1$OUTCOME == T)
X <- as.numeric(dat1$TREATMENT == T)

# Define other independent variables as one-hot coded data table

X <- 
  #mltools::one_hot() takes a data.table object instead of a tidy object
  mltools::one_hot(
    data.table::data.table(
      data %>%
        # Remember to remove treatment and outcome from X!
        dplyr::select(-OUTCOME, -TREATMENT)
    )
  )

set.seed(1714)

# Identify cases to be used for training and testing
# Training cases: random sample of 75% of observations in data

cases <- sample(seq_len(nrow(data)), round(nrow(data)*0.75))

# Split X, Y, and W separately by cases indices

train_X <- X[cases, ]
test_X <- X[-cases, ]

train_Y <- Y[cases]
test_Y <- Y[-cases]

train_W <- W[cases]
test_W <- W[-cases]

```

```{r Train CRF}

# Train Model 1

m1 <- grf::causal_forest(
  # Define Y, X, and W
  Y = train_Y,
  X = train_X,
  W = train_W,
  # Run more trees than default
  num.trees = 5000,
  # set seed internally
  seed = 1714
)

```

```{r CRF Treatment Effects}

priority.cate <- predict(
  object = m1,
  newdata = test_X,
  estimate.variance = T
)

test_X$preds <- priority.cate$predictions
test_X$Y <- test_Y
test_X$W <- test_W

```

```{r CRF Lift Scores}
# plot how sensitive each individual is to treatment

test_X %>%
  ggplot(aes(x = preds)) +
  geom_histogram() + 
  geom_vline(
    intercept = mean(train_Y[train_W == 1]) - mean(train_Y[train_W == 0]),
    col = 'black',
    lty = 'dashed') +
  labs(
    title = 'Distribution of Predicted Treatment Effects',
    subtitle = 'Average treatment effect in the training set shown with dashed reference line',
    y = 'Count',
    x = 'Predicted effect \np(OUTCOME|TREATMENT) - p(OUTCOME|CONTROL)') +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(text=element_text(family="serif"))

```

```{r CRF Plotting uncertainty}

# Plot predictions by their rank-order
ggplot(mapping = aes(
  x = data.table::frankv(preds$predictions, order = -1),
  y = preds$predictions))+
  geom_point() + 
  geom_errorbar(mapping = aes(
    ymin = preds$predictions + 1.96*sqrt(preds$variance.estimates),
    ymax = preds$predictions - 1.96*sqrt(preds$variance.estimates))) +
  geom_hline(yintercept = 0, col = 'black') +
  geom_hlin(
    yintercept = mean(train_Y[train_W == 1]) - mean(train_Y[train_W == 0]),
    lty = 'dashed', col = 'red') +
  labs(
    x = 'Rank',
    y = 'Estimated Treatment Effect',
    title = 'Predicted Individual Effects (Test Set)',
    subtitle = 'Estimate and 95% Confidence Interval',
    caption = 'Training set average treatment effect shown with dashed red line') +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(text=element_text(family="serif"))

```

```{r CRF RATE Test}

# Specify causal forest on the test (evaluation) set
m1_eval <- grf::causal_forest(
  Y = test_Y,
  X = test_X %>%
    #remove columns we previously added to the test_X
    dplyr::select(-preds, -Y, -W),
  W = test_W,
  num.trees = 5000,
  seed = 1714
)

# Compute RATE, using QINI

rate <- grf::rank_average_treatment_effect(
  forest = m1_eval,
  priorities = priority.cate$predictions,
  target = 'QINI',
)

# Compute a prioritization based on baseline risk
# Predicting the outcome as a function of independent variables, in the control condition
rf.risk <- regression_forest(X[cases, ][W[cases] == 0, ],
                             Y[cases][W[cases] == 0])
priority.risk <- predict(rf.risk, X[-cases, ])$predictions

# Test if two RATEs are equal
rate.diff <- rate_average_treatmen_effect(m1_eval,
                                          cbind(priority.cate$predictions,
                                                priority.risk))

rate.diff.tab <- data.frame(
  Target = rate.diff$target,
  Estimate = rate.diff$estimate,
  Lower = rate.diff$estimate - 1.96*rate.diff$std.err,
  Upper = rate.diff$estimate + 1.96*rate.diff$std.err
)

```

```{r Sources of Heterogeneity}

m1 %>%
  variable_importance() %>%
  as.data.frame() %>%
  mutate(Variable = colnames(m1$X.orig),
         V1 = round(V1, 3)) %>%
  arrange(desc(V1)) %>%
  rename(Important = V1) %>%
  slice(1:10) %>%
  knitr::kable(format = 'pipe', caption = 'Variable Important', label = 'grf_importance')

```
